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T. Grismayer^'QE. P. Alvefl 1 R. A. Fonseca, 1 ' 2 and L. O. SiW'Q 

1 GoLP/Instituto de Plasmas e Fusao Nuclear - Laboratorio Associado, 
Instituto Superior Tecnico, Lisbon, Portugal 
2 DCTI/ISCTE Instituto Universitdrio de Lisboa, 1649-026 Lisboa, Portugal 

(Dated: March 3, 2013) 

Abstract 

The generation of DC magnetic fields in unmagnetized plasmas with velocity shear is predicted 
for non relativistic and relativistic scenarios either due to thermal effects or due to the onset of 
the Kelvin-Helmholtz instability (KHI). A kinetic model describes the growth and the saturation 
of the DC field. The predictions of the theory are confirmed by multidimensional particle-in- 
cell simulations, demonstrating the formation of long lived magnetic fields (t ~ 100s along 
the full longitudinal extent of the shear layer, with transverse width on the electron length scale 
(\/7o c /u pe ), reaching magnitudes eB DC /m e ccj pe ~ /3 0v / 7o- 
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It is now recognized that kinetic plasma instabilities, operating on the electron time 
scale in unmagnetized plasmas, such as the Weibel instability [lj, can be of importance to 
explain the generation of magnetic fields relevant for shock formation in astrophysical and 
laboratory conditions and for particle acceleration and radiation emission in astrophysical 
scenarios [2j E] - Ab initio particle-in-cell simulations (PIC) have demonstrated the generation 
of subequipartition magnetic fields, i.e., the ratio between the energy density of the field and 
the kinetic energy density of the flow are close to 10 -3 — 10 -2 j4]-[6]. Previous studies have 
not considered the role of velocity shears in the dynamics of unmagnetized plasmas, which 
is known to lead to the onset of the collisionless Kelvin-Helmholtz instability (KHI). This 
scenario is of relevance from several perspectives: it has been proposed that the magnetic 
dynamo can operate over the seed magnetic field due to the collisionless KHI [3 [8] , it has 
been suggested that the KHI could play an important part in secondary magnetic island 
generation in magnetic reconnection sites [9j and the KHI remains a well known benchmark 
of magnetohydrodynamic simulations [TOj [11] . However the dynamics on the electron time- 
scale of the KHI is unexplored, and recent experiments have now addressed the collisional 
KHI [T2HT4] and in the near future, with the advent of more powerful lasers, experiments 
will be able to probe the collisionless KHI, in scenarios relevant for astrophysics. Motivated 
by recent particle-in-cell simulations that have demonstrated the onset of large scale DC 
magnetic fields in unmagnetized plasmas with velocity shear [15] , not captured by two-fluid 
theory [[5j [IS], we develop a theoretical model to predict the onset, growth and saturation 
of this magnetic field, as well as its spatial and long time behavior features. We compare the 
theoretical predictions with one (ID), two (2D) and three (3D) dimensional PIC simulations. 

We recall the 2D theoretical model of the unmagnetized KHI [[5], which is based 
on the relativistic fluid formalism of plasmas coupled with Maxwell's equations. Without 
loss of generality, we focus on symmetrically shearing flows (with velocities ±^ ^ along the 
y direction and with equal densities n ) with a tangential discontinuity in the x direction. 
The protons are considered free-streaming whereas the electron fluid quantities and fields are 
linearly perturbed, u = u e~ k± ^e l( ' k \\ y ~ u;t \ The unstable modes are stationary (Re (a;) = 0) 
surface waves and obey the following dispersion relation: 
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where V is the growth rate (Im(u;)) of the mode with wave number fey, uo pe = ^/47m e 2 /m e is 
the plasma frequency, k\ = k 2 +uo 2 e / (c 2 ^q)—uj 2 /c 2 and 70 = 1 / \J 1 — (3$ is the Lorentz factor 
of the shearing flows with /3 = vq/c. After a few e- folding times the system is dominated 
by the fastest growing mode with T max = ^/l/87^ 3 ^ 2 u pe and fc|| max = s/3/8^q 3 ^ 2 (jj pe /v . 

To ascertain these results and to fully explore the KHI, particle-in-cell (PIC) simulations 
were performed using OSIRIS [d fT8] . We simulate shearing slabs of cold (^ 3> ^/n where 
^t/i is the thermal velocity) unmagnetized electron-proton plasmas with a reduced mass ratio 
m p /m e = 100 (m e and m p are respectively the electron and the proton mass), and evolve 
it up to u) pi t = 100 (ujpi = ^m e /m p uj pe is the proton plasma frequency). First we present 
the simulation results of a shear flow with v = 0.2c. The shear flow initial condition is set 
by a velocity field with -\-v$e y in the middle-half of the simulation box, and a symmetric 
velocity field with —v$e y in the upper and lower quarters of the box. Initially, the system 
is charge and current neutral. Periodic boundary conditions are imposed in every direction. 
The simulation box dimensions are 500 x 200 (c/uj pe ) 2 , resolved with 20 cells per electron 
skin depth (c/oj pe ), and a number of 36 particles per cell per species are used. Space and 
time are respectively normalized to c/uj pe and l/uj pe - 

The growth and the wavenumber of the most unstable mode are confirmed by 2D simu- 
lations (Fig. [I] (al)). Interestingly, the growth of a DC (k\\ = 0) magnetic field mode is also 
observed (insets (a2) and (b2) of Fig. [T]), which is not predicted by the linear fluid theory 
(r(fc|| = 0) = 0) nor has it been previously identified in MHD simulations [HI [T9H2T] and 
only recently in kinetic simulations [15]. The growth of the DC magnetic field mode results 
from a current imbalance due to the mixing between the electron flows across the shear 
surface, while the proton flows remain almost unperturbed due to their inertia. The mixing 
arises due to the deformation of the electron interface between the two flows, which in the 
linearized fluid calculations is not accounted for and, in zeroth order, remains fixed. Alter- 
natively, we find that the physics describing the formation of a DC mode can be modeled in 
a ID reduced theory where an initial temperature drives the mixing effect. In order to un- 
derstand the formation of the DC field, we first consider the one-dimensional case. Initially 
all the fields are zero and we assume a plasma with a tangential shear flow in the x direc- 
tion and an initial temperature such that v t h <C ^o- It is clear that the thermal expansion 
of the electrons across the shear surface (ions are assumed to be cold and free streaming) 
leads to an imbalance of the current neutrality around the shear surface, forming a DC 

3 



y [C/G>pe] 

4 6 



10 



24 



25 



26 



al 



-v 




ct> pe t = 38 



-v 



B DC [m e c ft> pe /e] 
1.5xl0- : 



La2 


iii- 




, , , , i , , , , i , , , , i , , , , : 



4x10 

I- I 1 1 1 / I 1 1 1 I 1 1 1 I 1 - 



6xl0" 3 


-6xl0" 3 

io- 1 



3 

o 



PQ 




3 

o 



PQ 



1.5X10 1 





10 

120 

^30 

40 
50 



1 1 1 1 1 1 


1 | 1 1 1 | 1 


1 1 | 1 1 1 


1 cl 




-Vq 1 














1 co pe t = 1000 






- i i i 1 i i 


i 1 i i i 1 i 


1 1 1 1 1 1 - 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I . 

c2 




1.6x10" 



3 



1.6x10" 



FIG. 1. B z component of the magnetic field structure generated by the cold KHI for vq = 0.2c 
during (a) the linear regime, (b) near saturation and (c) at t = 1000 uj~} = 100 uj~1 . The insets 
on the right hand side represent the longitudinal average of the magnetic field, revealing the DC 
component. 



magnetic field in z direction. The initial corresponding electron distribution function reads 
f{x ) v X) v y , v Z) t = 0) = fo(v X i Vy — vo sgn(x), v z ). Due to the dimensionality of the problem, 
it is clear that E Z) B X) B y remain zero. The reduced set of equations is therefore Maxwell's 
equations coupled with the Vlasov equation d t F + v x d x F — e/m e (E + (v/c) x B z ).d$F = 
where F{x,v x ,v y ,t) = J dv z f{x,v x ,v y ,v z ,t). The formal solution of the Vlasov equation 
is F(x,v x ,v y ,t) = Fo(xo, ^xOj %))? where x , v x o and v y o denote the position and velocities 
of an electron at t — and F = J dv z ofo [22], [23]. At early times, if we assume that the 
induced fields are sufficiently small that we can neglect their effect on the change of momen- 
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turn of the electrons, the distribution function can be determined along the free streaming 
orbits of the electrons [24]. For the sake of simplicity, we separate the initial electronic 
distribution in two parts, F = F ~(xq < 0) + F + (x > 0). The electron currents read 
J^ y ~ —e f dVyVy J dv x F^(x — v x t,v x ,v y =p vq). For a Maxwellian distribution function, 
/mW = e~ v * /2v *h I y/2nv t h, we have F^(x ,v x0 ,v y0 =F v Q ) = n f M (v x0 )f M (v y0 =F v ) and we 
obtain J^ y ~ =Fen ^o J^ x / t dv x fM(v x ) — Te^ n erfc ^ - ^ . The total current is given by 
adding the unperturbed proton currents J v<y = en ^osng(x), from which the magnetic field 
can be integrated by neglecting the displacement current in Ampere's Law: 

B DC - e47m f3 V2v thx t 

where £ = \x\/y/2v t h x t- We thus verify that the DC magnetic field driven by the electron 
thermal expansion grows linearly with time. Its typical width is on the order of y2v t hxt and 
its peak intensity £> D c(x = 0) = 4v / ^en /3o^/^- This derivation is valid as long as the 
orbits of the electrons do not diverge much from the free streaming orbits, i.e., as long as 
the electric and magnetic fields that develop self-consistently do not affect the free motion 
of the particles. However, the electrons will eventually feel the induced magnetic field which 
tends to push more electrons across the shear via the v x £? D c force. Consequently, the rate 
at which electrons cross the shear increases, which, in turn, enhances the growth rate of the 
magnetic field. 

In order to verify our analytical calculations and to further investigate the phase where 
electrons deviate from their free streaming orbits, we carried out ID simulations of the 
x direction of the 2D simulations. The Debye length is resolved in the ID simulations 
(Ax = Xd) and we used 1000 particles per cell. Fig. [2] (al-3) shows the time evolution 
of the xp x phase space and the magnetic field for v t h = 0.016c and v = 0.2c. At earlier 
times uj pe t = 9 (Fig. [2] (al)) an excellent agreement between the model and the simulation 
is observed. The model breaks down approximately when an electron initially with v x = 
(around the shear) acquires a velocity change on the order of Vthx, which corresponds to a 
strong distortion of the Maxwellian distribution around the shear. Fig. [2] (a2) shows this 
effect at uj pe t = 17. The model underestimates the magnitude of the magnetic field and one 
can clearly observe the distortion of the distribution function in the field region. As the 
magnetic field grows, the Larmor radius (tl) of the electrons crossing the shear interface 
decreases. When the minimum r^min (associated to the peak of £> D c) becomes smaller than 
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the characteristic width of the magnetic field foe? the bulk of the electrons becomes trapped 



trapping prevents the electron bulk expansion across the shear (that drives the growth of 
the magnetic field), saturating the magnetic field. An estimate of the saturation can be 
obtained by equating r^min ~ foe- From Eq. ([2]), it is possible to write the magnetic field as 
Bvc(x,t) = ^7reno(3ow(x,t), where w(0,t) should be interpreted as the characteristic width 
of the field. With Z DC ~ w(Q,t), r L5min = m^o7o/em e S DC (0, t), we find that Z D c ~ c-s/l^/^pe 
giving the saturation level of the magnetic field as eB^/m e cuj pe ~ Poy/To- This scaling has 
been verified for ID simulations (Fig. [3]). 

In the absence of an initial temperature, an alternative mechanism is needed to drive 
the electron mixing across the shear surface that in turn generates the DC field. This 
mechanism is the cold fluid KHI that operates in 2D and 3D geometries. In fact, in the 
warm shear flow scenario, both the cold fluid KHI and the electron thermal expansion can 
contribute to the generation of the DC field. This happens when the typical length of the 
DC field due to the thermal expansion (Zdc) after a few e- foldings of the cold fluid KHI 
(T K Hi- g rowth = ^e-foidings/r ma x, where n e _ fo i dings is on the order of 10) is on the order of the 
relativistic electron skin depth, i.e., ^/i^KHi-growth ~ s/%c/uj pe . Therefore, the cold fluid 
KHI dominates the electron mixing in the limit ^Ikm-growth <C /^ pe . 

For a two dimensional cold plasma undergoing the KHI, the electron distribution function 
can be written as f(x,y,v x ,v y ,v z ,t) = n S(v x - v xf i(x,y,t))6(v y - v yf i(x,y,t))6(v z ) where 
v xfi>Vyfi correspond to the velocity field solutions of the fluid theory. In this case, the self- 
generated KHI fields play the role of an effective temperature that transports the electrons 
across the shear surface, while the protons remain unperturbed, inducing a DC component 
in the current density, and hence in the fields. We then have to solve the evolution of 
the distribution function and show that the current density J yi averaged over a wavelength 
A = 27r/fc||, has a non zero DC part. We follow the same approach as before and calculate 
the average distribution function defined as: 



To obtain analytical results we will assume that the linearly perturbed fluid quantities are 
purely monochromatic, which is equivalent to assume that after a few e-foldings, the mode 
corresponding to fen = fc|| max dominates with a growth rate of Y — Y max . We then write 



by the magnetic field structure. This is illustrated in Fig. [2](a3) at uj pe t = 55. The magnetic 
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Vyfi ~ Vq(x) and v x fi = v x fi sm(k\\y)e~ k± \ x \ +rt , where v x fi, the amplitude of the velocity 
perturbations at t — 0, is associated to the small thermal fluctuations (small enough to 
ensure that the thermal expansion is negligible over T K Hi- g rowth)- Inserting v x fi,v y fi into 
Eq. ([3]), we obtain 

F{*,v x ,t) = ^==, (4) 

where £(x,v x ,t) = v x /vmax(%it) with v max (x,t) = v x fie~ k± ^ +Tt . We observe that the devel- 
opment of 2D cold KHI reveals close similarities with the ID hot model previously described. 
In the 2D KHI, averaging the distribution in the direction of the flow shows that the per- 
turbation gives rise to a spread in v x that may be interpreted as an effective temperature. 
The spread in v x decays exponentially away from the shear and grows exponentially with 
time. The mean velocity is zero and the effective temperature associated to this distribution 
function is defined as V£ff(x,t) = (l/^o) / dv x v 2 x F(x,v xl i) = v% nax /2. One can then expect 
a similar physical picture as in the hot shear scenario and, as a result, the emergence of 
DC components in the fields which are induced by the development of the unstable KH 
perturbations. The evolution of the phase space in Fig. [2] illustrates the similarity between 
the warm ID (insets al-3) and cold 2D (insets bl-3) scenarios. 

The challenge in this scenario is to determine how such a distribution function expands 
across the shear surface due to the complexity of the orbits in the fields structure (multi- 
dimensional fields with discontinuities at x — 0). In the region where the electron mixing 
occurs, we assume electron orbits given by x ~ x + (v x o/T)e rt and v x ~ v x $e vt where x 
and v x0 are the position and velocity of a particle at the time t when the instability begins. 



J max 



J e, y ( x ^) — Tev / dv x F(x,v x ,t) (5) 
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where xT G [— ^ max , ^ max ] an d ^ max (^) = ^max(^ = 0, t) that represents the maximum velocity 
of a particle that was originally in the vicinity of the shear. The limits of the integral 
([5]) represent the deformation of the boundary between the two flows on a characteristic 
distance of ^ max /r as the instability develops. In the fluid theory, the boundary remains 
fixed, precluding the development of the DC mode. We then find the total current density 



by summing the proton contribution and integrate to obtain the induced DC magnetic field: 



B B c(x ^ 0,t) = =F47ren A) 

The peak of the DC magnetic field is located at x = where the expression above reduces to 
£>dc(0, t) = 8e^o^o^max(^)/( 7I_ r) and thus grows at the same rate as the KHI fields. One can 
verify in Fig. [2] (bl-b2) that Eq. shows reasonable agreement with the 2D simulations. 
This derivation neglects the DC Lorentz force on the electron trajectories, which makes this 
model valid as long as induced DC field remain small compared to the fluid fields associated 
to the mode k\\ max . The peak of £>dc field is proportional to ^ max (£). We therefore conclude 
that the induced DC magnetic field is always on the same order of the fluid fields (Fig. [j] 
(a) and (b)) and thus its consequences to KHI development cannot be neglected. As the 
DC field evolves, electrons start to get trapped and we expect a level of saturation similar 
to the one obtained in the ID model which is verified by the simulations. The comparisons 
between the saturation level of the ID, 2D and 3D simulations are shown in Fig. |3j verifying 
the Poy/7o scaling. Interestingly, the DC magnetic field remains stable beyond the electron 
time scale and persists up to 100s uj~^ as shown in Fig. [I] (c). Eventually the protons will 
drift away from the shear surface due to the magnetic pressure, broadening the DC magnetic 
field structure and lowering its magnitude. 

In conclusion, we have presented an analytical description of the formation of a DC 
magnetic field in a shear flow scenario, which is in good agreement with ID, 2D and 3D 
PIC simulations. We have shown that the DC magnetic field can arise due to an electron 
thermal expansion across the shear, in a warm shear scenario; we also extended this picture 
to the cold shear scenario where the development of the cold fluid KHI perturbations induce 
an effective temperature that drives the electron mixing across the shear. The DC magnetic 
field saturates on the electron time scale and persists up to proton time-scales, reaching 
maximum magnitudes of eB^Qjmcijj^ ~ 1.4/3 y7o with thicknesses of a few sf%c/uj pe) and 
thus is dynamically relevant for the evolution of the KHI on ion time scales. 
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Simulations were performed at the 1ST cluster (Lisbon, Portugal), and the Jugene super- 
computer (Germany). 
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FIG. 2. Evolution of the electron phasespace. Left: ID warm shear flow with vq = 0.2 c and 
^th — 0.016. Right: 2D cold shear flow with vq = 0.2c. The blue (red) color represents the electrons 
with a negative (positive) drift velocity vo. The self-consistent DC magnetic field is represented 
by the solid curve, whereas the dashed curve represents the magnetic field given by the theoretical 
model. 
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FIG. 3. Magnitude of the DC magnetic field peak at saturation as a function of /?o\/7o- The red, 
blue and black markers represent the results of ID, 2D and 3D PIC simulations, respectively. The 
error bars are associated to the fluctuations of the peak value in the saturation stage. The lines 
represent best fit curves to the simulation results, demonstrating the good agreement with Eq. ??. 
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